nanopyx.core.analysis.cross_correlation_elastic

View Source

  1import numpy as np
  2from math import sqrt
  3from skimage.filters import gaussian
  4from scipy.interpolate import interp2d
  5
  6from .ccm import calculate_ccm_from_ref
  7from .estimate_shift import GetMaxOptimizer
  8from ..transform.blocks import assemble_frame_from_blocks
  9
 10
 11def calculate_translation_mask(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel", algorithm="weight"):
 12    
 13    if algorithm == "weight":
 14        return calculate_translation_mask_vector_weight(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method=method)
 15    elif algorithm == "field":
 16        return calculate_translation_mask_vector_field(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method=method)
 17    else:
 18        print("Not a valid algorithm option! Select either 'weight' or 'field'")
 19
 20def calculate_translation_mask_vector_weight(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel"):
 21
 22    width = img_slice.shape[1]
 23    height = img_slice.shape[0]
 24
 25    assert width == img_ref.shape[1]
 26    assert height == img_ref.shape[0]
 27
 28    block_width = int(width / blocks_per_axis)
 29    block_height = int(height / blocks_per_axis)
 30
 31    flow_arrows = []
 32    blocks_stack = []
 33
 34    for y_i in range(blocks_per_axis):
 35        for x_i in range(blocks_per_axis):
 36            x_start = x_i * block_width
 37            y_start = y_i * block_height
 38
 39            slice_crop = img_slice[y_start:y_start+block_height, x_start:x_start+block_width]
 40            ref_crop = img_ref[y_start:y_start+block_height, x_start:x_start+block_width]
 41            slice_ccm = np.array(calculate_ccm_from_ref(np.array([slice_crop]).astype(np.float32), np.array(ref_crop).astype(np.float32))[0])
 42
 43            if max_shift > 0 and max_shift*2+1 < slice_ccm.shape[0] and max_shift*2+1 < slice_ccm.shape[1]:
 44                ccm_x_start = int(slice_ccm.shape[1]/2 - max_shift)
 45                ccm_y_start = int(slice_ccm.shape[0]/2 - max_shift)
 46                slice_ccm = slice_ccm[ccm_y_start:ccm_y_start+(max_shift*2), ccm_x_start:ccm_x_start+(max_shift*2)]
 47
 48            if method == "subpixel":
 49                optimizer = GetMaxOptimizer(slice_ccm)
 50                max_coords = optimizer.get_max()
 51                ccm_max_value = -optimizer.get_interpolated_px_value(max_coords)
 52            else:
 53                max_coords = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
 54                ccm_max_value = slice_ccm[max_coords[0], max_coords[1]]
 55
 56            ccm_width = slice_ccm.shape[1]
 57            ccm_height = slice_ccm.shape[0]
 58            blocks_stack.append(slice_ccm)
 59
 60            if ccm_max_value >= min_similarity:
 61                vector_x = (ccm_width/2.0 - max_coords[1] - 1)
 62                vector_y = (ccm_height/2.0 - max_coords[0] - 1)
 63                flow_arrows.append([x_start + block_width/2.0, y_start + block_height/2.0, vector_x, vector_y])
 64
 65    if len(flow_arrows) == 0:
 66        print("Couldn't find any correlation between frames... try reducing the 'Min Similarity' parameter")
 67        return None
 68
 69    translation_matrix = np.zeros((height, width*2))
 70    translation_matrix_x = np.zeros((height, width))
 71    translation_matrix_y = np.zeros((height, width))
 72
 73    max_distance = sqrt(width * width + height * height)
 74
 75    for j in range(height):
 76        for i in range(width):
 77            # iterate over vectors
 78            dx, dy, w_sum = 0, 0, 0
 79
 80            if len(flow_arrows) == 1:
 81                dx = flow_arrows[0][2]
 82                dy = flow_arrows[0][3]
 83
 84            else:
 85                distances = []
 86                all_distances = 0
 87                for arrow in flow_arrows:
 88                    d = sqrt(pow(arrow[0] - i, 2) + pow(arrow[1] - j, 2)) + 1
 89                    distances.append(d)
 90                    all_distances += pow(((max_distance - d) / (max_distance * d)), 2)
 91
 92                for idx, arrow in enumerate(flow_arrows):
 93                    d = distances[idx]
 94                    first_term = pow(((max_distance - d) / (max_distance * d)), 2)
 95                    second_term = all_distances
 96
 97                    weight = first_term / second_term
 98                    dx += arrow[2] * weight
 99                    dy += arrow[3] * weight
100                    w_sum += weight
101
102                dx = dx / w_sum
103                dy = dy / w_sum
104
105            translation_matrix_x[j, i] = dx
106            translation_matrix_y[j, i] = dy
107
108    if blocks_per_axis > 1:
109        translation_matrix_x = gaussian(translation_matrix_x, sigma=max(block_width, block_height/2.0))
110        translation_matrix_y = gaussian(translation_matrix_y, sigma=max(block_width, block_height/2.0))
111
112    translation_matrix[:, :width] += translation_matrix_x
113    translation_matrix[:, width:] += translation_matrix_y
114
115    blocks = assemble_frame_from_blocks(np.array(blocks_stack), blocks_per_axis, blocks_per_axis)
116
117    return translation_matrix, blocks
118
119def calculate_translation_mask_vector_field(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel"):
120    """
121    Function used to calculate a translation mask between 2 different images.
122    Method based on dividing both images in smaller blocks and calculate cross correlation matrix between corresponding
123    blocks. From the ccm, the shift between the two images is calculated for each block and a translation matrix is
124    using the shifts in the center position of each block and then interpolating the remaining translation mask.
125    :param img_slice: numpy array with shape (y, x); image to be used for translation mask calculation
126    :param img_ref: numpy array with shape (y, x); image to be used as reference for translation mask calculation
127    :param max_shift: int; maximum shift accepted between each corresponding block, in pixels.
128    :param blocks_per_axis: int; number of blocks for both axis
129    :param min_similarity: float; minimum similarity (cross correlation value after normalization) between corresponding
130    blocks.
131    :param method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve
132    subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.
133    :return: numpy array with shape (y, x) where width is equal to two times the width of the original image.
134    [:, :width/2] corresponds to the translation mask for x and [:, width/2:] corresponds to the translation mask for y.
135    """
136
137    width = img_slice.shape[1]
138    height = img_slice.shape[0]
139
140    assert width == img_ref.shape[1]
141    assert height == img_ref.shape[0]
142
143    block_width = int(width / blocks_per_axis)
144    block_height = int(height / blocks_per_axis)
145
146    blocks_stack = []
147
148    y_translation = []
149    x_translation = []
150
151    for y_i in range(blocks_per_axis):
152        for x_i in range(blocks_per_axis):
153            x_start = x_i * block_width
154            y_start = y_i * block_height
155
156            slice_crop = img_slice[y_start:y_start+block_height, x_start:x_start+block_width]
157            ref_crop = img_ref[y_start:y_start+block_height, x_start:x_start+block_width]
158            slice_ccm = np.array(calculate_ccm_from_ref(np.array([slice_crop]).astype(np.float32),
159                                                        np.array(ref_crop).astype(np.float32))[0])
160
161            ccm_x_start = 0
162            ccm_y_start = 0
163
164            if max_shift > 0 and max_shift*2+1 < slice_ccm.shape[0] and max_shift*2+1 < slice_ccm.shape[1]:
165                ccm_x_start = int(slice_ccm.shape[1]/2 - max_shift)
166                ccm_y_start = int(slice_ccm.shape[0]/2 - max_shift)
167                slice_ccm = slice_ccm[ccm_y_start:ccm_y_start+(max_shift*2), ccm_x_start:ccm_x_start+(max_shift*2)]
168
169            if method == "subpixel":
170                optimizer = GetMaxOptimizer(slice_ccm)
171                max_coords = optimizer.get_max()
172                ccm_max_value = -optimizer.get_interpolated_px_value(max_coords)
173            else:
174                max_coords = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
175                ccm_max_value = slice_ccm[max_coords[0], max_coords[1]]
176
177            blocks_stack.append(slice_ccm)
178
179            if ccm_max_value >= min_similarity:
180                shift_x, shift_y = get_shift_from_ccm_slice(slice_ccm, method=method)
181                y_translation.append([y_start + max_coords[0] + ccm_y_start,
182                                      x_start + max_coords[1] + ccm_x_start,
183                                      shift_y - 0.5])
184                x_translation.append([y_start + max_coords[0] + ccm_y_start,
185                                      x_start + max_coords[1] + ccm_x_start,
186                                      shift_x - 0.5])
187
188    y_translation = np.array(y_translation)
189    x_translation = np.array(x_translation)
190    y_interp = interp2d(y_translation[:, 0], y_translation[:, 1], y_translation[:, 2])
191    x_interp = interp2d(x_translation[:, 0], x_translation[:, 1], x_translation[:, 2])
192
193    translation_matrix = np.zeros((height, width*2))
194    translation_matrix_x = np.zeros((height, width))
195    translation_matrix_y = np.zeros((height, width))
196
197    for j in range(translation_matrix_x.shape[0]):
198        for i in range(translation_matrix_x.shape[1]):
199            translation_matrix_x[j, i] = x_interp(j, i)
200            translation_matrix_y[j, i] = y_interp(j, i)
201
202    translation_matrix[:, :width] += translation_matrix_x
203    translation_matrix[:, width:] += translation_matrix_y
204
205    blocks = assemble_frame_from_blocks(np.array(blocks_stack), blocks_per_axis, blocks_per_axis)
206
207    return translation_matrix, blocks
208
209def get_shift_from_ccm_slice(slice_ccm, method="subpixel"):
210    """
211    Function used to calculate the shift corresponding to the maximum correlation between two images.
212    :param slice_ccm: numpy array with shape (y, x);
213    :param method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve
214    subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.
215    :return: tuple of floats; values corresponding to shift_x and shift_y, in this order.
216    """
217
218    w = slice_ccm.shape[1]
219    h = slice_ccm.shape[0]
220
221    radius_x = w / 2.0
222    radius_y = h / 2.0
223
224    if method == "subpixel":
225        optimizer = GetMaxOptimizer(slice_ccm)
226        shift_y, shift_x = optimizer.get_max()
227    elif method == "Max":
228        shift_y, shift_x = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
229
230    shift_x = radius_x - shift_x - 0.5
231    shift_y = radius_y - shift_y - 0.5
232
233    return (shift_x, shift_y)

def calculate_translation_mask( img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method='subpixel', algorithm='weight'): View Source

12def calculate_translation_mask(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel", algorithm="weight"):
13    
14    if algorithm == "weight":
15        return calculate_translation_mask_vector_weight(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method=method)
16    elif algorithm == "field":
17        return calculate_translation_mask_vector_field(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method=method)
18    else:
19        print("Not a valid algorithm option! Select either 'weight' or 'field'")

def calculate_translation_mask_vector_weight( img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method='subpixel'): View Source

 21def calculate_translation_mask_vector_weight(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel"):
 22
 23    width = img_slice.shape[1]
 24    height = img_slice.shape[0]
 25
 26    assert width == img_ref.shape[1]
 27    assert height == img_ref.shape[0]
 28
 29    block_width = int(width / blocks_per_axis)
 30    block_height = int(height / blocks_per_axis)
 31
 32    flow_arrows = []
 33    blocks_stack = []
 34
 35    for y_i in range(blocks_per_axis):
 36        for x_i in range(blocks_per_axis):
 37            x_start = x_i * block_width
 38            y_start = y_i * block_height
 39
 40            slice_crop = img_slice[y_start:y_start+block_height, x_start:x_start+block_width]
 41            ref_crop = img_ref[y_start:y_start+block_height, x_start:x_start+block_width]
 42            slice_ccm = np.array(calculate_ccm_from_ref(np.array([slice_crop]).astype(np.float32), np.array(ref_crop).astype(np.float32))[0])
 43
 44            if max_shift > 0 and max_shift*2+1 < slice_ccm.shape[0] and max_shift*2+1 < slice_ccm.shape[1]:
 45                ccm_x_start = int(slice_ccm.shape[1]/2 - max_shift)
 46                ccm_y_start = int(slice_ccm.shape[0]/2 - max_shift)
 47                slice_ccm = slice_ccm[ccm_y_start:ccm_y_start+(max_shift*2), ccm_x_start:ccm_x_start+(max_shift*2)]
 48
 49            if method == "subpixel":
 50                optimizer = GetMaxOptimizer(slice_ccm)
 51                max_coords = optimizer.get_max()
 52                ccm_max_value = -optimizer.get_interpolated_px_value(max_coords)
 53            else:
 54                max_coords = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
 55                ccm_max_value = slice_ccm[max_coords[0], max_coords[1]]
 56
 57            ccm_width = slice_ccm.shape[1]
 58            ccm_height = slice_ccm.shape[0]
 59            blocks_stack.append(slice_ccm)
 60
 61            if ccm_max_value >= min_similarity:
 62                vector_x = (ccm_width/2.0 - max_coords[1] - 1)
 63                vector_y = (ccm_height/2.0 - max_coords[0] - 1)
 64                flow_arrows.append([x_start + block_width/2.0, y_start + block_height/2.0, vector_x, vector_y])
 65
 66    if len(flow_arrows) == 0:
 67        print("Couldn't find any correlation between frames... try reducing the 'Min Similarity' parameter")
 68        return None
 69
 70    translation_matrix = np.zeros((height, width*2))
 71    translation_matrix_x = np.zeros((height, width))
 72    translation_matrix_y = np.zeros((height, width))
 73
 74    max_distance = sqrt(width * width + height * height)
 75
 76    for j in range(height):
 77        for i in range(width):
 78            # iterate over vectors
 79            dx, dy, w_sum = 0, 0, 0
 80
 81            if len(flow_arrows) == 1:
 82                dx = flow_arrows[0][2]
 83                dy = flow_arrows[0][3]
 84
 85            else:
 86                distances = []
 87                all_distances = 0
 88                for arrow in flow_arrows:
 89                    d = sqrt(pow(arrow[0] - i, 2) + pow(arrow[1] - j, 2)) + 1
 90                    distances.append(d)
 91                    all_distances += pow(((max_distance - d) / (max_distance * d)), 2)
 92
 93                for idx, arrow in enumerate(flow_arrows):
 94                    d = distances[idx]
 95                    first_term = pow(((max_distance - d) / (max_distance * d)), 2)
 96                    second_term = all_distances
 97
 98                    weight = first_term / second_term
 99                    dx += arrow[2] * weight
100                    dy += arrow[3] * weight
101                    w_sum += weight
102
103                dx = dx / w_sum
104                dy = dy / w_sum
105
106            translation_matrix_x[j, i] = dx
107            translation_matrix_y[j, i] = dy
108
109    if blocks_per_axis > 1:
110        translation_matrix_x = gaussian(translation_matrix_x, sigma=max(block_width, block_height/2.0))
111        translation_matrix_y = gaussian(translation_matrix_y, sigma=max(block_width, block_height/2.0))
112
113    translation_matrix[:, :width] += translation_matrix_x
114    translation_matrix[:, width:] += translation_matrix_y
115
116    blocks = assemble_frame_from_blocks(np.array(blocks_stack), blocks_per_axis, blocks_per_axis)
117
118    return translation_matrix, blocks

def calculate_translation_mask_vector_field( img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method='subpixel'): View Source

120def calculate_translation_mask_vector_field(img_slice, img_ref, max_shift, blocks_per_axis, min_similarity, method="subpixel"):
121    """
122    Function used to calculate a translation mask between 2 different images.
123    Method based on dividing both images in smaller blocks and calculate cross correlation matrix between corresponding
124    blocks. From the ccm, the shift between the two images is calculated for each block and a translation matrix is
125    using the shifts in the center position of each block and then interpolating the remaining translation mask.
126    :param img_slice: numpy array with shape (y, x); image to be used for translation mask calculation
127    :param img_ref: numpy array with shape (y, x); image to be used as reference for translation mask calculation
128    :param max_shift: int; maximum shift accepted between each corresponding block, in pixels.
129    :param blocks_per_axis: int; number of blocks for both axis
130    :param min_similarity: float; minimum similarity (cross correlation value after normalization) between corresponding
131    blocks.
132    :param method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve
133    subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.
134    :return: numpy array with shape (y, x) where width is equal to two times the width of the original image.
135    [:, :width/2] corresponds to the translation mask for x and [:, width/2:] corresponds to the translation mask for y.
136    """
137
138    width = img_slice.shape[1]
139    height = img_slice.shape[0]
140
141    assert width == img_ref.shape[1]
142    assert height == img_ref.shape[0]
143
144    block_width = int(width / blocks_per_axis)
145    block_height = int(height / blocks_per_axis)
146
147    blocks_stack = []
148
149    y_translation = []
150    x_translation = []
151
152    for y_i in range(blocks_per_axis):
153        for x_i in range(blocks_per_axis):
154            x_start = x_i * block_width
155            y_start = y_i * block_height
156
157            slice_crop = img_slice[y_start:y_start+block_height, x_start:x_start+block_width]
158            ref_crop = img_ref[y_start:y_start+block_height, x_start:x_start+block_width]
159            slice_ccm = np.array(calculate_ccm_from_ref(np.array([slice_crop]).astype(np.float32),
160                                                        np.array(ref_crop).astype(np.float32))[0])
161
162            ccm_x_start = 0
163            ccm_y_start = 0
164
165            if max_shift > 0 and max_shift*2+1 < slice_ccm.shape[0] and max_shift*2+1 < slice_ccm.shape[1]:
166                ccm_x_start = int(slice_ccm.shape[1]/2 - max_shift)
167                ccm_y_start = int(slice_ccm.shape[0]/2 - max_shift)
168                slice_ccm = slice_ccm[ccm_y_start:ccm_y_start+(max_shift*2), ccm_x_start:ccm_x_start+(max_shift*2)]
169
170            if method == "subpixel":
171                optimizer = GetMaxOptimizer(slice_ccm)
172                max_coords = optimizer.get_max()
173                ccm_max_value = -optimizer.get_interpolated_px_value(max_coords)
174            else:
175                max_coords = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
176                ccm_max_value = slice_ccm[max_coords[0], max_coords[1]]
177
178            blocks_stack.append(slice_ccm)
179
180            if ccm_max_value >= min_similarity:
181                shift_x, shift_y = get_shift_from_ccm_slice(slice_ccm, method=method)
182                y_translation.append([y_start + max_coords[0] + ccm_y_start,
183                                      x_start + max_coords[1] + ccm_x_start,
184                                      shift_y - 0.5])
185                x_translation.append([y_start + max_coords[0] + ccm_y_start,
186                                      x_start + max_coords[1] + ccm_x_start,
187                                      shift_x - 0.5])
188
189    y_translation = np.array(y_translation)
190    x_translation = np.array(x_translation)
191    y_interp = interp2d(y_translation[:, 0], y_translation[:, 1], y_translation[:, 2])
192    x_interp = interp2d(x_translation[:, 0], x_translation[:, 1], x_translation[:, 2])
193
194    translation_matrix = np.zeros((height, width*2))
195    translation_matrix_x = np.zeros((height, width))
196    translation_matrix_y = np.zeros((height, width))
197
198    for j in range(translation_matrix_x.shape[0]):
199        for i in range(translation_matrix_x.shape[1]):
200            translation_matrix_x[j, i] = x_interp(j, i)
201            translation_matrix_y[j, i] = y_interp(j, i)
202
203    translation_matrix[:, :width] += translation_matrix_x
204    translation_matrix[:, width:] += translation_matrix_y
205
206    blocks = assemble_frame_from_blocks(np.array(blocks_stack), blocks_per_axis, blocks_per_axis)
207
208    return translation_matrix, blocks

Function used to calculate a translation mask between 2 different images. Method based on dividing both images in smaller blocks and calculate cross correlation matrix between corresponding blocks. From the ccm, the shift between the two images is calculated for each block and a translation matrix is using the shifts in the center position of each block and then interpolating the remaining translation mask.

Parameters

img_slice: numpy array with shape (y, x); image to be used for translation mask calculation
img_ref: numpy array with shape (y, x); image to be used as reference for translation mask calculation
max_shift: int; maximum shift accepted between each corresponding block, in pixels.
blocks_per_axis: int; number of blocks for both axis
min_similarity: float; minimum similarity (cross correlation value after normalization) between corresponding blocks.
method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.

Returns

numpy array with shape (y, x) where width is equal to two times the width of the original image. [:, :width/2] corresponds to the translation mask for x and [:, width/2:] corresponds to the translation mask for y.

def get_shift_from_ccm_slice(slice_ccm, method='subpixel'): View Source

210def get_shift_from_ccm_slice(slice_ccm, method="subpixel"):
211    """
212    Function used to calculate the shift corresponding to the maximum correlation between two images.
213    :param slice_ccm: numpy array with shape (y, x);
214    :param method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve
215    subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.
216    :return: tuple of floats; values corresponding to shift_x and shift_y, in this order.
217    """
218
219    w = slice_ccm.shape[1]
220    h = slice_ccm.shape[0]
221
222    radius_x = w / 2.0
223    radius_y = h / 2.0
224
225    if method == "subpixel":
226        optimizer = GetMaxOptimizer(slice_ccm)
227        shift_y, shift_x = optimizer.get_max()
228    elif method == "Max":
229        shift_y, shift_x = np.unravel_index(slice_ccm.argmax(), slice_ccm.shape)
230
231    shift_x = radius_x - shift_x - 0.5
232    shift_y = radius_y - shift_y - 0.5
233
234    return (shift_x, shift_y)

Function used to calculate the shift corresponding to the maximum correlation between two images.

Parameters

slice_ccm: numpy array with shape (y, x);
method: str, either "subpixel" or "max"; defaults to "subpixel"; subpixel uses a minimizer to achieve subpixel precision in shift calculation. max simply takes the coordinates corresponding to the max value of the ccm.

Returns

tuple of floats; values corresponding to shift_x and shift_y, in this order.